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Abstract 

We present a new Monte Carlo algorithm applying a history-driven mechanism for the calculation 
of the density of states for classical statistical models. With the new method, detailed balance is 
naturally satisfied in limit and the estimated density of state converges to the exact value. The 
new method could be easily evolved into the multicanonical method to achieve high accuracy. 

PACS numbers: 64.60.Cn,05.50.+q,02.70.Lq 



Among the various algorithms proposed to obtain the density of states {g{E)) through 
simulation [H-E], one effective strategy is to apply a history-driven mechanism to the 
Metropolis sampling Monte Carlo (MC) simulation [JHO]. The idea is to use the infor- 
mation earned from the history of the simulation (e.g., the histogram H{E) accumulated) 
as a feedback to determine the transition probability of the future moves and hence the gen- 
eration of future states, so that the simulation is driven to sample the entire (or targeted) 
energy range efficiently and provides an estimation of the density of states. In general, a MC 
trial move from energy Ei to E2 could be accepted according to the Metropolis-Hastings 
transition probability: 

p{E2^E,)=mm[l,R{H)], (1) 

where R(H) is a certain function of the histogram and the core of the history-driven method. 
Starting with a random initial guess, the history-driven simulation goes like: 

HjE,) _ giEME2 ^ E,) _ gjE,) 

H{E,) g{E2)p{E, ^ E2) g{E2) ^ ^' ^ ^ 

and the density of states that is in priori unknown could be estimated as: 

Certain correlation between the states generated in near sequence is necessary for a history- 
driven mechanism to take effect, if, otherwise, the generation of a new state is completely 
independent to the history of the simulation, there wont be any history-driven mechanism. 

If R{H) is fixed during the simulation, Eq. (2) describes the well known multicanonical 
simulation. A simple history-driven mechanism is the multicanonical recursion [7j. For 
example, the entropic sampling jl] adopts recursive periods of multicanonical simulations: 

Starting with a random initial guess, such as Ri{H) = 1, the transition probability is to be 
updated at the end of each period using the histogram accumulated from that period: 

To make the method practical, the histogram is initialized to be Hi{E) = /iq < 1 at the 
beginning of each period 
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In practice, each period of the entropic simulation must be long enough so that a useful 
working estimation of the density of states could be obtained for a final long run to achieve 
high accuracy. However, the long periods also slow the updating of Rn{H) and reduce the 
efficiency of the simulation in sampling the targeted energy range, hindering the application 
of this method to a relatively large system. A more timely updated transition probability 
could drive the simulation to go through the energy range faster. Obviously, the limiting 
case of doing this is to update R{H) right after every single MC move. With this single-move 
multicanonical recursion, Eq. (5) turns into: 

R^H) = l^j^j = /^(^^)-^(^i), (6) 

where H{E) is the histogram maintained throughout the simulation and / > 1 is a modi- 
fication factor. (Note that Hi{E) = ho initially and Hi{Ei) = ho + 1 if Ei is visited after 
the only trail move of the z*'^ period, so |g| = (^) ^ Similarly, |^ = ^ if ^2 
is visited.) This is nothing but the transition probability used in the well-known Wang- 
Landau (WL) algorithm [S], which is actually a limiting case of multicanonical recursion [S]. 
With a properly chosen modification factor, usually f = e = 2.71828, the WL algorithm is 
very efficient in sampling throughout the energy range, conquering the probably dramatic 
energy landscape of the system. As noted before, e.g., [12|, the estimation of density of 
states |[g = = fH{E,)-H{E2) (note that the WL algorithm results into a fiat 

histogram) dose not converge to the exact value with a fixed /. High accuracy could also be 
achieved for the WL algorithm by doing "global updates" to R{H) recursively (in addition 
to the "local updates" that are taken after every single move): 

Rn{H) = i?„_i(/7)/„^"(^^)-^"(^^) (7) 

where /„ is carefully reduced every time the histogram Hn{E) is fiat enough. To avoid 
saturation of errors, we use the ln(/) = 1/t scheme [12] for the WL algorithm in comparing 
with the new method in the rest of this paper. 

Not confining to the multicanonical recursion, we may go back to Eq. (2) and ask that 
whether other form of R{H) could be found as an effective history-driven mechanism to 
achieve high efficiency and accuracy. We find that R{H) = is a fairly good choice, 

where a is properly chosen power index. Applying the new method, the simulation goes like: 

HjE,) _giE,) ( H{E^) Y 

H{E2) g{E2) [h{E,)J ' ^ ^ 
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Detailed balance is naturally satisfied in limit and the histogram converges to: 



(9) 



H{E2) \gmj 



In other words, the density of states could be estimated as: 



) 



a+l 



(10) 



which converges to the exact value without any extra efforts of doing global updates. 

According Eq. (9), a sufficiently large power index a should be used to make the histogram 
relatively flat and achieve efficiently sampling throughout the targeted energy range. A 
simple choice of a is to make exp(a) to be about the total number of configurations of the 
system as that ensures a flatness of ^ e. Tested for a 16 x 16 2D Ising model with a 

total number of configurations of 2^^^ ~ e^'''^, we find that any a > 200 works almost as well 
in driving the simulation to sample throughout the entire energy space using ~ 5.0±1.0 x 10^ 
flips, while a = 100 doubles the simulation effort with ~ 9.0±2.0x 10^ flips. As a comparison, 
the WL algorithm also takes ~ 5.0 ± 1.0 x 10^ flips with ln{f) > 0.1 and ~ 8.0 ± 2.0 x 10^ 
flips with ln{f) = 0.05. It is worth to point out that, since the new method is as easily 
realized as the WL algorithm (see below), measurement in MC trial moves (flips) is almost 
equivalent to that in CPU time in comparing these two methods. So the new method is as 
efficient as the WL algorithm in sampling the energy space not only in MC trial moves but 
also in real CPU time. 

Furthermore, we test the efficiency of the new method for a traveling salesman problem 
(TSP) consisting of 100 cities. For 100 points ("cities") randomly distributed within a 
1x1 square, the distances [E in this problem) of the shortest and the longest cyclical path 
connecting the 100 cities found after different simulation efforts are given in Tab. |T} Trail 
moves proposed by [Tl] are used and the path length is binned into bins of width 0.01 for 
the maintenance of H{E) in the simulation. Fixed value of / = and a = 100000 is 
used for the WL and the new method, respectively. We repeat the tests for several different 
randomly generated distributions of the cities and find similar results as shown in Tab. |T} 
Again, the new method is as efficient as the WL method. 

With the muticanonical method being a special case (a = 0) of the new method, the 
latter shares many characters of the former. First of all, the errors of the estimated density 
of states do not saturate even without doing global updates (i.e., with a fixed power index 
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a). So the new method could be useful for the estimation of the density of states in the 
cases where global updates are prohibited or not so feasible. For example, both the WL 
simulation [H] and the new method simulation could be used to find out the energy space and 
the possible energy entries of a system, which, in most not theoretically known 

as is the case of Ising model. For this purpose, the parameter / or a shouldn't be reduced 
during the simulation, and the new method is superior to the WL method by providing 
an estimation of the density of states that converges to the exact value. Moreover, the 
histogram resulted from the new method is not necessary to be flat. Actually, with a proper 
guess of g{E), the histogram could be roughly designed for various purposes. 

The new method is also closely related to the WL method. Given the well known result 

HiE) ^ 

of the harmonic series that for HiE) — )■ oo, 'S^ — = In H(E) + c, where c is the Euler's 

^1 ^ 

constant , we have 

Comparing to Eq. (6), we can see that only a small modification to the WL algorithm is 
needed to implement the new method, i.e., instead of doing In g{E) In g{E) +ln/, we do 
lng{E) ^ lng{E) + j^j^ for the local updates. The histogram H{E) must also be maintained 
throughout the simulation, but with negligible calculating effort. Also, we see that the the 
power index a here is an equivalent of the modification factor In / in WL method. Actually, 
for the 2D Ising model, we find that the fluctuation of the histogram is proportional to :^7=^, 
consistent with the finding in [13] that the fluctuation of the histogram resulted from WL 
simulation is proportional to i — . From Equ. (10), we can predict that the error in the 



estimated density of states is proportional to ^^=pj = \/a + 1, which is exactly what we see 
in Fig [l] in long run. 

If a better estimation of the density of states is available, a smaller power index a could 
be used to achieve high accuracy with reduced simulation effort. For example, using the 

n I M 

initial guess g{E) = g{M) = Cn ^ for the 2D Ising model (where n is the total number 
of spins and M the magnetization), a could an order of magnitude smaller than that needs 
for the initial guess g{E) = 1. So we can reduce a recursively when better estimation of the 
density of states is available after a period of simulation and eventually turns the simulation 
into a multicanonical one. We here propose a scheme of global updates of R{H) that results 
into both high efficiency and accuracy. After the first period of simulation that samples the 



targeted energy range, we update the transition probabihty Ri{H) = ( ^[^^j ) globally into 

where b C [0, 1] is a modification factor. The histogram H{E) is maintained throughout 
the simulation and whose value at the end of the n*'* period is Hn{E). In general, after 
the entire energy space is sampled by a period of simulation, the following global update is 
taken, 

Rn+i{H) = J„MH) ff^l , (13) 



where 

J„+i(i7) = UH) [^;^^) (14) 

and Ji{H) = 1. For this purpose, a separated histogram Hp{E) is maintained for each 
period and we apply a criterion Hp{E) > Hmin for ending a period of simulation and doing 
global updates. Testing for the Ising model, we find that a modification factor of 6 = 0.5 
works robustly with an Hmin — I x I, where / is the size of the Ising model. Seen from 
Eq. (11), by doing the global updates, the new method simulation seamlessly evolvesinto 
a muticanonical one as 6" — )■ and the second term on the right hand side vanishes when 
n goes large. Actually, we terminate the global update and explicitly turn the simulation 
into a muticanonical one when ab"' < 0.1 in our simulations. ' We can see from Fig. [TJthat 
after applying the global updates, the new method is as accurate as the 1/t WL method. In 
Fig. |2| we also show the specific heat C{t) calculated from the g{E) for the 1024 x 1024 Ising 
model, which are obtained from simulations of the same length but applying the new method 
or the WL method. Again, the new method is as accurate as the 1/t WL method. We also 
point out that the global update scheme we proposed here is even more easily realized than 
the 1/t scheme in the WL method as no calculation of t is needed. 

Comparing to the 1/t WL method, the only disadvantage we find for the new method is 
that it does not try to maintain a flat histogram self-adaptively. The histogram may get not 
so flat and reduce the efficiency of simulation in long run. But we find this problem can be 
easily solved by using proper values of a, b, Hmin- The guide line is to use a larger b or Hmin 
for larger a to ensure the fiatness of the histogram, but may in the cost of slightly slower 
convergency, as seen from Fig. |2} 
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The new method is extremely general and many techniques developed based on the 
muticanonical method and/or the WL method could be applied to the new method. 

We thank D. P. Landau, Shan-Ho Tsia, Meng Meng, and Chenggang Zhou for helpful 
discussion. This work is partly supported by NSFC grant No. 11143008. 
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TABLE I: Minimum and maximum path length found by the new method (NM) and the WL 
method simulations for a TSP that consists of 100 cities. The average path length and the error 
bars are calculated from 20 independent runs. See text for more details. 





Min 


Max 


MC moves 


NM WL 


NM WL 


10^ 


23.7 ±0.43 23.7 ±0.40 


74.5 ±0.22 74.4 ±0.17 


10^ 


12.7 ±0.23 12.6 ±0.15 


77.38 ±0.08 77.39 ±0.06 


10^ 


8.13 ±0.20 8.13 ±0.06 


77.81 ±0.01 77.81 ±0.01 
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1/t Wang-Landau method: fi=e 




t (MC moves) 



FIG. 1: Comparison of the accuracy of the new method and the 1/t WL method for g{E) of 
the 16 X 16 2D Ising model. No sign of saturation of errors is found for the new method even 
without doing global updates {b = 1.0). Applying the global updates mentioned in the text, the 
new method is as accurate as the 1/t WL method. Each curve is averaged from 100 independent 
runs. The error in the density of states is calculated as a[lgg{E)] = \j ^ ^'^ (-^) L ^ where 
ge{E) is the exact value and the summation is taken over all the energy entries of total number n. 
The normalization g{{)) = 3e(0) = In 2 is used for the calculation. 
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FIG. 2: Specific heat for the 256 x 256 2D Ising model. The relative errors (e(C) = | ^ |, where 
Ce is the exact value) are shown in the inset. Both the new method simulation and the WL method 
simulation make use of a single run of 5.12 x 10^^ MC moves. 
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